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We explore the numerical stability properties of an evolution system suggested by Alekseenko and 
Arnold. We examine its behavior on a set of standardized testbeds, and we evolve a single black 
hole with different gauges. Based on a comparison with two other evolution systems with well- 
known properties, we discuss some of the strengths and limitations of such simple tests in predicting 
numerical stability in general. 
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I. INTRODUCTION 

In recent years the quest for finding a numerically sta- 
ble formulation of the Einstein evolution equations has 
become more and more intense, see e.g. for reviews. 
Effort has been put into finding first order symmetric 
hyperbolic formulations of the evolution equations, since 
the properties of such systems can be analyzed mathe- 
matically. However, even if an evolution system is sym- 
metric hyperbolic there is no guarantee that its numerical 
implementation is stable when evolving a highly dynamic 
black hole spacetime. The Baumgarte-Shapiro-Shibata- 
Nakamura (BSSN) [il3 system is an example of a system 
that is not first order symmetric hyperbolic but has nice 
stability properties. Thus, current mathematical analysis 
is not sufficient to explore the properties of an evolution 
system. The system must be implemented and tested nu- 
merically before we can draw definite conclusions about 
its viability for a particular physical application. 

In 5] , Alekseenko and Arnold ( AA) suggest a first or- 
der formulation of the evolution equations that is sym- 
metric hyperbolic when considering only a subset of the 
variables. In particular, the metric itself and some of 
its first spatial derivatives are not considered part of the 
evolution system and are treated as given functions when 
showing hyperbolicity. It is argued that this is sensible 
since the metric is derivable from an ordinary differential 
equation. A distinguishing feature of the AA system is 
that a minimal number of first derivatives of the met- 
ric are introduced as independent variables. The system 
has only 20 unknowns and no parameters that have to 
be fixed for hyperbolicity, so it is relatively simple for a 
symmetric hyperbolic system. 

We have chosen to implement the AA system numeri- 
cally and to compare its numerical properties with those 
of the Arnowitt-Deser-Misner (ADM) H, Q and BSSN 
systems. The ADM system is known to be unstable 
in many situations, but we include it here since for fi- 
nite time intervals of evolutions it has been used rather 
successfully in practice , e . g. i n 3D black hole simula- 
tions H 11 [lOl In m 111 Til m. The BSSN system 
can be obtained from a trace-conformal decomposition 



of the ADM equations, and with suitable techniques it 
is very stable for black hole evolutions. For example, 
the first stable evolution for more than 100000 A/ of a 
single Schwarzschild black hole in a (3-f l)-dimensional 
code without adapted coordinates was obtained with a 
modified BSSN system Both the ADM and BSSN 
systems are not first order. However, there exist first 
order versions of the BSSN system which are symmet- 
ric hyperbolic 0, 0, llH |]3 > if the densitized lapse and 
shift are considered given functions. Straightforward first 
order forms of ADM are only weakly hyperbolic, which 
implies certain numerical instabilities (see e.g. |2C|). 

The second order version of the BSSN system shares 
with the AA system the property that a subsystem of it is 
indeed symmetric hyperbolic. A crucial step in the con- 
struction of BSSN is the introduction of the contracted 
Christoffel symbol of the conformal metric, F*, as an in- 
dependent variable. The evolution equation for the ex- 
trinsic curvature then contains derivatives of the metric 
only in the form of a Laplace operator, so that the met- 
ric obeys a wave equation if F' is considered a prescribed 
variable, ignoring its own evolution equation. This par- 
tial hyperbolicity of the BSSN system may well be a cru- 
cial ingredient in its success as evolution system. Hence 
the question arises whether symmetric hyperbolicity in 
a subsystem implies a numerical advantage for the AA 
system. 

There is a variety of numerical tests that one can per- 
form on an evolution system. A complex and important 
issue is that seemingly minor changes in the test condi- 
tions and implementation of the system can lead to very 
different conclusions about stability. In [23|, an impor- 
tant step is taken towards creating a set of tests that 
can serve as a standardized benchmark for stability. Our 
discussion of the AA system in comparison to ADM and 
BSSN can also be viewed as a contribution to the ongoing 
development of such benchmarks. On the one hand, we 
report on the performance of our particular implemen- 
tation of ADM and BSSN, where a large body of prior 
work allows us to judge how representative the current 
benchmark is. On the other hand, we apply these tests 
to a new evolution system about which nothing is known 
so far from numerical experiments, and the question is. 
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for example, whether one can predict the usefulness of 
the AA system for black hole evolutions. We therefore 
also discuss results for the evolution of a single black hole 
that go beyond the current set of benchmarks in ^| . 

The paper is organized as follows. In Scc.m we write 
out the AA system explicitly, in notation that is more 
familiar to numerical relativists, and we show that the 
complete system is not symmetric hyperbolic. We intro- 
duce a time-independent conformal rescaling of the AA 
system that we will need later for black hole evolutions. 
In Sec. mil we perform three tests from the recently pro- 
posed test suite for numerical evolutions, the robust 
stability test, the gauge stability test, and the linear wave 



stability test. Finally, we evolve a single black hole space 
time in Sec. II VI since this is the physics example that we 
are most interested in. We conclude with a discussion in 
Sec. El 



II. THE AA SYSTEM 

In 0, the evolution system is given in compact nota- 
tion. Here we write it out in the form in which wc have 
implemented it: 



/yfc = TTlKidik,] - 9ok,t + {{gp3,q- gpq,i) gik- {gpi,q- gpq,i) gjk) g^'') , (1) 



w 



ijk 
7) 



2\/2 

^(/5',5p. + /3^.<?P,), (2) 
K = .g'^/Cy, (3) 

\ {'^9pq.jgZ - gpq,.gZ + 2g„,.5^', - 2g^,., + " g"' g'l) - (4) 

gij,sgpq,rg^'^ g + gpq,rgsi,jg^'^ g + gpq.rgsj^ig^'^ g ) i 

S = \ (2clj + 2^] ~ gp,^qg'^%9r^g'"' + gpq^^f^grrg"" - gpr,qgZ9rof' + 5P9.«/!.5.,5P^+ (5) 

gpq,rg''jgsig^'' + gpq,rg'''tgsjg'"' - gpq,Tg'''jgs^g'^'' - gpq^rg'^'igsjg^'' - gpj,qg'^'',sgrtg'^''- 
gpqjg'^^grzg^' - gpuqg'^^gng^' - gpq,ig'''^,sgr3g''' + 2gpj,q5^^s5r»5''' + ^gpi^^gP'^^gr^g'''- 
'^gr3.sgpq,rg^'' g"^" + 4:gij,sgpq,rgP''g'"') , 

- (2acfj ~ 2a J - a^pgij^qg'^'^ + a^pg^^jg^"^ + a^pgqj^^gP'' + 2aKKij+ (6) 
2(3P^Kp, + 2(3P^Kp, - A.agP'^Kp.Kq,) , 

■'^ = -l{gpq,r{2g'^W-gZr + 9st.,ugP'g''''g'')) (7) 

i (2c^j - a (g.3,sgpq..r (- i9''9'l + g^'g^ + g^J (c'^ KpqKP'^))) , (8) 

igrj.pgtk - grt.pgjk) g'^'' + /3 ^ {gik,p - gqr.pgikg'''') + /3^, {-9jk,p + 9qr,pgjkg'^'')) , (9) 

^ {'^^fjk + \/2 (m,fcj - Wjka + Wpj.qgikg^'^ - Wpq^gitg'"^ - Wpi,qg3kg^'' + Wpq^igjkg'"'- (10) 

" {9p3.q9^''Ktk + gpq.jgP'^K.k + gpi^qg^'^K^k - gpq,ig^'^Kjk + gpj,qgikKP'^ - gpq.jgikK^'^- 
gpi.qgjkKP'' + gpq,igjkKP'') + gpj.qg^'^w^k - gpqjg^'^Wik - gpi.qg'"'wjk + gpq.ig^'^Wjk- 
gpj.^gikw^'^ + gpq.jgikw^'^ + gpi.qgjkw'^'^ - gpq.igjkw^'')) , 

Cijk = - (4cJ]j, + V2 {2Kajgik - 2Ka^igjk - 2ag'"jqg3kKpi + 2agP'jggikKpj - ag'P''^gikKpq+ (11) 

<^9^'!i9jkKpq + 2a,pgjkg'"^Kqi - 2a,pgikg^'' Kqj - agpq^rgjkg^'^g'''' K^i + agpq.rgikg^'' g'''' Ksj)) , 
c^o5y = -2aKij +2wij, (12) 

d^K,3 = \ (2i?y + V2a (g^"^ (/p., -f /„,) + /p.,, + - g''%fps3gq.g'' + /'./.,p5,./^~ (13) 

g^irfpsigqjg'''' + g^^fsipgqjg'''')) , 

dofijk = Cijk + {-aKikj + aKjk,t - ajKik + a^iKjk) ■ (14) 
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The operator do = dt — P^dp is defined in terms of ordi- 
nary partial derivatives for all variables. 

The new variable fijk is anti-symmetric in i and j and 
also satisfies the cyclic property /yfc -I- fjki + fkij = 0- 
This means that fijk represents 8 independent variables. 
The total number of 20 variables is thus smaller than in 
the case of most first order formulations (which usually 
involve 30 or more variables), but larger than for the 
BSSN system, which has 15 independent variables. 

In addition, the equations in the AA system are more 
complex than in the BSSN system, so that the AA evolu- 
tion system takes roughly twice as long to run the same 
number of iterations, even though we have made an ef- 
fort to implement the AA system in a numerically op- 
timal way. The simulations were carried out with the 
BAM code, which is a rewritten version of the code used 
in . Part of BAM is a Mathematica script to convert 
tensor equations to C code, and the AA equations were 
generated this way using Mathematica and MathTensor. 



A. Hyperbolicity of the AA system 

Looking at Eqs. l(T|l- H14|) and assuming that a and /3* 
are prescribed given functions, we see that the system 
has the form 

dtu + A\u)u,i + v{u)+w{m''^{u)u^^Q{u)uj) =0. (15) 
Here 




(16) 



is the solution vector with spatial indices suppressed, and 





A'{u) = ( r'{u) a'{u) c'{u) 
^s'{u) c'{u) b\u)^ 

is a matrix which contains a symmetric submatrix 

^a*(u) c'(u)^ 



S\u) 



d{u) V{u) 



(17) 



(18) 



In addition, Eq. (|15|l contains the two vector valued func- 
tions V and w, where the argument of w depends on 



(19) 



and its transpose U"^, and also on another matrix Q{u), 
which is of the simple form 




^q{u) 0^ 
Qiu) = I I , 
Oy 




so that the argument of w in fact only depends on squares 
of first derivatives of gij . 

From Eq. H15|) it is immediately apparent that the sys- 
tem is of first order form, as no second derivatives of u 
appear. Nevertheless the system differs from the stan- 
dard first order form by the term w. This means that 
the standard theorems about well-poscdness (see e.g. Q) 
are not applicable, and that we cannot easily compute 
characteristic speeds and modes of the system. However, 
if we linearize the system around any background , 
all terms of the form w drop out. This can be seen as 
follows. Assume that 



(21) 



where is a small perturbation to the background . 
Then 

w (m^-' {u)u^^Q{u)u 

= w{m'^iu)g^,qiu)g^^) + 

w {m^={u){g^,q{u)g^^ + 5,i9(«).9,f )) + O {{g'^f) 
= w {m^={u)g^,q{u)g^) + d\u^)g^, + O {{g'^f) ,(22) 

and hence Eq. H15|) becomes 

dtu^ + A^u^ -f v{u^) + O {{u^ f) + w{u^) = 0, (23) 

which is now of standard first order form, but with a 
modified matrix 





a*(u^) c'(u^) 



(24) 



(20) 



Looking at Eq. H23|) and H24|) , it is immediately clear that 
the system is not symmetric hyperbolic. In fact, since 
the matrix S^{u^) has several zero eigenvalues, the full 
system in general is unlikely to have a complete set of 
eigenvectors and thus to be strongly hyperbolic. 

The exception is linearizion around flat space where P 
and P in Eq. H24|l vanish, so that the system H23|l be- 
comes symmetric hyperbolic in this special case. Also, 
since S'^{u) is symmetric, the subsystem of Kij and fijk 
with gij considered as a prescribed variable is symmetric 
hyperbolic, which holds true even for the non-linearized 
system (|15|l . Alekseenko and Arnold Q emphasize this 
last point and stress that the evolution equation for gij 
is (as usual) just an ordinary differential equation. Yet 
in numerical evolutions the latter property may be unim- 
portant, since gij has to be evolved along with the other 
variables and cannot be prescribed. In addition, it is not 
clear if the fact that the evolution system (|15|l contains a 
symmetric hyperbolic subsystem has any bearing on the 
stability properties of the full system. 
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B. Conformal version of the A A system 



III. STABILITY TESTS 



In Sec. IIVI we report on numerical evolutions of black 
hole puncture data 0, . In order to numerically 

evolve puncture data without excision the AA system has 
to be modified. The reason is that e.g. for two punctures 
the metric has the form 



with conformal factor 



mi 
2^1 



9tj 

1712 

' 2^ 



■u, 



(25) 



(26) 



where ri and r2 are distances from the punctures, mi and 
TO2 are the bare puncture masses, and u is finite. From 
Eq. (|25|l it is apparent that the physical metric diverges 
like 



1 

9^3 ~ ^ 



(27) 



at each puncture, so that finite differencing calculations 
across or near any puncture will fail. The same problem 
also occurs in the variable 



(28) 



For this reason we do not evolve the variables g^j and 
fijk directly, rather we rescale gij and fijt by the time 
independent conformal factor 



mi mo 

V' = 1 + — + — 

2ri 2r2 
such that the rescaled quantities 



9^j 

fij k 



(29) 



(30) 
(31) 



become finite at the puncture. Furthermore, we also in- 
troduce the rescaled extrinsic curvature 



(32) 



K,, 



(33) 



Since t/; is an a priori prescribed time independent func- 
tion, the principal part of the system of evolution equa- 
tions (fTT)- (|14|l remains unchanged if we use the rescalings 
in Eqs. H3Q|I . H31|l . and H32|l to express the system in terms 
of the new variables gij, Kij, and fijk- In addition, all 
terms involving g.^, Kij, fijk, and their derivatives are 
finite in the new system so that finite differencing can 
be used without trouble. There are, however, additional 
terms containing spatial derivatives of "0, which cannot 
be computed using finite differencing. Yet, since ip is & 
known function we can use analytic expressions for its 
derivatives. Furthermore, all such derivative terms also 
contain appropriate powers of V' which make them finite. 



We have used some of the stability testbeds suggested 
in [23, also known as the "apples with apples" tests, 
to explore the properties of the AA system. We also 
show test results for the ADM and BSSN systems for 
comparison. 



A. Robust stability test 

The purpose of the robust stability test is to determine 
how an evolution system will handle random errors that 
are bound to occur at machine-precision level. Random 
constraint violating initial data in the linearized regime 
is used to simulate this machine error. 

All apples with apples test are run on a full-3D grid, 
but in this case with periodic boundary conditions in all 
directions. The 3D domain is a 3-torus, and here we 
only use 3 grid-points in the y- and ^-directions. The 
parameters are: 

• Simulation domain: x G [—0.5, +0.5] 

• Number of grid points in each direction: = 50p, 
Uy — Uz — 3, with p ~ 1,2,4 

• Courant factor = 0.5 

• Iterations: lOOOOOp, output every lOOp iterations 
(corresponding to 1 crossing time) 

• Gauge: harmonic, i.e. dta — —a^trK, /3' = 
The initial data is given by 



-ij , 



a — 

= 



(34) 
(35) 
(36) 
(37) 



in order to remove divergences in Kij , which in the case where e. 
of punctures with spin behaves as 



a random number with a probabil- 



ity distribution, which is uniform in the interval 
(-10"i"/p2,+10"i7p^)- After a smaU number of 
timesteps, the random noise in gij will have propagated 
into all other evolved quantities, except for which 
will remain identically zero for all time. Hence our ini- 
tial data differ slightly form the ones in who add 
random noise to all quantities which need initialization. 

There is one obvious problem with the robust stability 
test. The initial data are random, and we must carefully 
check that the random number generator we use does 
not introduce errors because the numbers are not totally 
random. We use C's "random" function on Linux Red- 
hat 7.3, a pseudo-random number generator based on a 
non-linear additive feedback algorithm that avoids the 
short-comings of some of the older implementations of 
the "rand" function. As a seed for the random number 
generator we use the time where the subroutine is called 
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plus the clock cycle of the processor on which it is called. 
Since our code is parallelized and each processor uses its 
own seed, the actual random numbers on the grid depend 
on the number of processors. This could be avoided by 
additional coding, however, the qualitative result of the 
robust stability test is not expected to depend on the ac- 
tual random numbers. We have run the test several times 
on different numbers of processors. We have also tried 
increasing the size of the domain, both in the x-direction 
and the y- and z-directions (which results in a different 
number of random numbers being generated). We find 
that although this does change the results, many features 
are robust and do not change with these variations. We 
will only comment on these robust features of the test. 

ADM and BSSN results are similar to those reported 
in I23. ADM grows exponentially (Fig. P), while BSSN 
is stable (Fig. EJ. Note that ADM clearly shows the sig- 
nature of a grid mode, since when plotting versus the 
number of iterations almost identical exponential growth 
is obtained for the three different resolutions. 

In 24] , the ADM system is run with a Courant factor of 
only 0.03 with the result that it is stable for much longer, 
up to 200 crossing times. We have lowered the Courant 
factor in the ADM run to 0.25 and found that we do not 
see exponential growth in the Hamiltonian constraint. 
When plotting tiK we see some growth, similar to Fig. 2 
in [231, but we do not encounter a blowup (we ran the 
coarsest resolution to 10000 crossing times, and did not 
encounter exponential growth). The results are shown 
in Fig. |2] We also tried with Courant factor 0.03 as 
in |24| . and this did not show any exponential blowup 
(the run was stopped at 600 crossing times). Preliminary 
experiments indicate that the transition between stable 
and unstable Courant factors for ADM lies between 0.4 
and 0.5. This deserves further investigation. 

The results for AA are shown in Fig. 01 AA is quali- 
tatively stable for the 1000 crossing times suggested for 
this test, but oscillations in the maximum and minimum 
of the Hamiltonian occur. There are two types of oscil- 
lations, one with a long wavelength and one with a short 
wavelength. We can see that the short wavelength oscil- 
lations dampen out. The long wavelength of the medium 
and high resolution runs seems to be the same, while 
the wavelength in the low resolution run is about half of 
the wavelength in the other runs. As pointed out earlier, 
these features are robust when changing the domain-size, 
but the amplitude of the oscillations are affected by this 
change. 

Given the presence of these slow oscillations, we ran 
another set of tests where we evolve for 10000 crossing 
times. Fig. [S] shows the results. We can see that at late 
times, the constraint violation for all 3 resolutions are 
growing, and this suggests that the AA system is stable 
for random noise initial data for a long time, but it will 
eventually become unstable. We also ran the test of the 
BSSN system to 10000 crossing times but found no insta- 
bilities or any indication of longterm growing oscillations 
in this case. 



We tried lowering the Courant factor for the AA runs. 
We ran the coarsest resolution to 10000 crossing times 
with a Courant factor of 0.25 and found that the growth 
in the Hamiltonian constraint still happens, but it hap- 
pens later than in the run with the Courant factor of 0.5. 
Lowering the Courant factor for the AA system does not 
change the general features of the oscillations. 



B. Gauge wave stability test 

In this test we look at the ability of the evolution sys- 
tems to handle gauge dynamics. This is done by consid- 
ering flat Minkowski space in a slicing where the 3-metric 
gij is time dependent. The gauge wave is then given by 

/27r(a;-<)\ , , 

= 1 + Asm^- ^ , (38) 



9yy ~ 9zz — Ij 
9xy — 9xz — 9yz ~ O7 



K — 



ttA 



COS 



1 + A sm ( — — '- 



K, 



otherwise 



(39) 
(40) 

,(41) 
(42) 



a = ,11 + Asm { ^^^^ ), (43) 



(44) 



Here d is the size of the domain in the i-direction and A 
is the amplitude of the wave. Since this wave propagates 
along the x-axis and all derivatives are zero in the y- and 
z-directions, the problem is essentially one dimensional. 

These are the parameters we use in our gauge stability 
tests: 

• Simulation domain: x G [—0.5, +0.5] 

• Points: rix = 50p, riy — — 3, p — 1,2,4 

• Courant factor — 0.25 

• Iterations: 200000p, output every 200p iterations 
(corresponding to 1 crossing time) 

• ^ = 0.1 and A = 0.01 

• Gauge: harmonic, i.e. dta — —a^tiK, (3^—0 

Since we know the analytical solution at all times, we 
can compare our numerical results to the analytical re- 
sults, see also 20]. As well as testing if the system has 
exponentially growing modes, we can check the conver- 
gence of the numerical result to the analytical solution 
with increasing resolution. 

Fig.Elshows that for a gauge wave amplitude A = 0.1 
the AA system converges as expected for a finite time 
interval, but there is exponential growth and the run 
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crashes after about 100 crossing times. For A = 0.01 
(Fig. [7|) it takes longer for this non-convergence to show, 
about 1000 crossing times, but it is stiU there. If the 
AA system together with the given gauge equations was 
symmetric hyperbohc, and if one had a stable discretiza- 
tion scheme, then the result would be convergent. Con- 
versely, assuming that ICN is an appropriate discretiza- 
tion scheme, we would conclude that the AA system 
in harmonic gauge is not symmetric hyperbolic, which 
agrees with our result in Sec. Ill Al 

The ADM system shows no growth of the constraints 
at all, and the constraint violation remains at machine 
precision. This is probably because for a 1-dimensional 
gauge wave the ADM system simplifies such that no con- 
straints violating modes are possible. 

The results for the BSSN system are shown in Fig. |S1 
The BSSN system becomes unstable and non-convergent 
rather quickly (on the order of 100 crossing times). This 
is somewhat surprising since the BSSN system has been 
very successful in single black hole evolutions and neutron 
star evolutions [23ll25l |. so we expect it to be able to han- 
dle gauge dynamics. However, part of the robustness of 
BSSN can be attributed to the fact that constraint violat- 
ing modes have a finite speed and can propagate 
out of the grid for example for radiative boundary con- 
ditions. However, on our grid with periodic boundaries, 
when a constraint violating mode appears it will stay on 
the grid, which is probably the reason that the BSSN 
system fails this test. Thus, this test by itself cannot 
determine whether a system can handle gauge dynamics, 
but must be accompanied by other tests without peri- 
odic boundaries. In particular, there is no contradiction 
to the observed stability of BSSN for single black holes 
with radiative boundaries. 



C. Linear wave stability test 



given by 



An ( 2'Kx 
At: ( lux 



(47) 

(48) 
(49) 



d 

Kij = 0, otherwise. 
These are the parameters of our run: 

• Simulation domain: x G [—0.5, +0.5] 

• Points: — 50p, Uy — Uz — 3, p — 1,2,4 

• Courant factor — 0.25 

• Iterations: 200000p, output every 200p iterations 
(corresponding to 1 crossing time) 

• A^ 10-8 

• Gauge: geodesic, i.e. a = 1,/?' = 0. 

Figs. lUITUl show the results for the ADM, AA, and 
BSSN systems, respectively. Fig. |21 is a comparison of 
the numerical wave to the analytical wave at 100 cross- 
ing times. We see that the AA numerical wave travels 
slightly faster than the ADM and BSSN numerical waves. 
Fig. El shows the L2-norm of the difference between the 
numerical and analytical linear waves as a function of 
time, and again the wave in the AA system travels faster 
than in the other systems. The Hamiltonian constraint 
has a value of about 10"^ at the end of the run for the AA 
system, so we are still well within the linear regime. It is 
also worth noting that in the ADM system the constraint 
violation is constant throughout the evolution, while the 
constraints grow slightly for the other two systems (but 
the maximum violation is of the order 10^^ for the entire 
run, for all resolutions). 



In this section we investigate the ability of the evolu- 
tion systems to propagate the amplitude and phase of 
a traveling gravitational wave. The amplitude of this 
wave is sufficiently small so that we are in the linear 
regime. This test reveals effects of numerical dissipation 
and other sources of inaccuracy in the evolution algo- 
rithm. 

The initial 3-metric and extrinsic curvature are ob- 
tained from 

ds^ = -df + dx^ + (1 + 6) dy^ + (1 - 6) (^^^ (45) 
where 

b^Asini—^ ^J, (46) 

d is the size of the propagation domain, and A is the 
amplitude of the wave. The extrinsic curvature tensor is 



IV. NUMERICAL TESTS INVOLVING A 
SINGLE SCHWARZSCHILD BLACK HOLE 

In this Section we present numerical results for the evo- 
lution of a single Schwarzschild black hole in two different 
gauges. 

A. Geodesic slicing 

The initial data for this test is a Schwarzschild black 
hole in isotropic coordinates. We use geodesic slicing 
to evolve the initial data. In geodesic slicing, the coor- 
dinate lines correspond to freely falling observers, and 
the resulting slicing of Schwarzschild can be expressed 
analytically in terms of Novikov coordinates, e.g. [28l| . 
Chapter 2.7.2, and 29], which we transform to isotropic 
coordinates for direct comparison with the numerical re- 
sults. To this end, we numerically solve the following 
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implicit equation for the Schwarzschild area radial coor- 
dinate R — R{t, i?inax), 




(50) 



where R = i?max is the position at r = of a freely falling 
observer starting at rest, and M is the mass of the black 
hole. To transform from the i?max radial coordinate to 
the isotropic radial coordinate r we use 



{M + 2r) 
4r 



(51) 



Then the rr-component of the metric in isotropic coor- 
dinates is given by 



grr = i^{r) 



4 / OR 



\dR„ 



(52) 



where 



dR 



1 R 



dRn 



2R„ 



R 



1 arccos 




(53) 



and R is computed from Eq. 15U|) . Here '0(r) = 1 + l^- 
Note that there is a typo in 29], Eq. (16). 

Analytically, the metric becomes infinite at time 



tr 



ttM, 



(54) 



which is the time for an observer that starts from rest at 
the horizon to reach the Schwarzschild singularity. This 
results in a crash in the numerical computations. We run 
this test on a 3D grid in the so-called cartoon mode [33| , 
because the spherical symmetry of the problem means 
that we can do a computation using information only 
on a line passing trough the center of the black hole. 
The computational domain is z S [0, 80M],x = 0, y = 0. 
We use second order accurate finite differencing together 
with an iterative Crank-Nicholson scheme with a Courant 
factor of 0.25 for evolution. Since we never run longer 
than tcrash = ttA/, the outer boundary at 80M has no 
effect on the black hole located at the origin. 
In order to compute the order of convergence. 



a = log2 



fih — f: 



2h 



f2h — fh 



(55) 



we run the test at the resolutions h = O.OIM, 2h = 
0.02M, and 4/i = 0.04M. We compute a on each grid- 
point present in the coarsest run, and fh is the value of 
the quantity under consideration (here the metric com- 
ponent gzz) for resolution h. Our grid-points are not in 



the same places for all three resolutions, because we al- 
ways stagger the puncture between two grid-points. This 
means that we must interpolate to get the values at all 
the coarse grid-points, for which we use 4th order La- 
grange interpolation. Fig. II II shows a for the AA system 
at time t = and time t — 3.0Af in the region close to 
the black hole where the metric deviates the most from 
the flat conformal metric. We see that we have second 
order convergence close to the black hole in both cases. 
The spikes in the curve at later times are due to the fact 
that the curves of gzz for the 3 resolutions cross. 



B. 1+log lapse, gamma driver shift 

We have implemented l-flog lapse and gamma driver 
shift in our code. This gauge choice makes the BSSN sys- 
tem stable for more than lOOOAf for certain single black 
hole runs |2^. While the geometric motivation of this 
gauge choice, namely singularity avoidance and reduc- 
tion of slice-stretching, should be valid for any evolution 
system, it is unclear whether the AA system will be as 
stable as BSSN or unstable in this test case. In par- 
ticular, note that the gamma driver shift is designed to 
control the evolution of one of the BSSN variables, F*. 
We have not implemented proper outer boundary condi- 
tions for the AA system, so in our simulations we have 
waves coming in from the outer boundaries. We use the 
same parameters as for the geodesic runs, i.e. the domain 
in cartoon mode is z G [Q,80M],x — 0,y = 0, the res- 
olution is ft, = O.OlAf, and we use a Courant factor of 
0.25. 

Our AA run crashes at around 20 Af. Very steep gra- 
dients appear near the black hole and they eventually 
kill the run. The reason for this seems to be that the 
1+log lapse depends on the trace of the physical extrin- 
sic curvature. Since we evolve the conformal extrinsic 
curvature, finite difference errors are enlarged by a factor 
of ip'^ when computing the physical extrinsic curvature, 
and these account for the differences in the lapse evolu- 
tion between BSSN and AA runs. When using the AA 
system, the lapse drops very fast, leaving a sharp gradi- 
ent between the frozen region and the region where the 
fields can evolve. The gradients in the physical variables 
that kill the run appear where the frozen region meets 
the dynamic region. In the BSSN system this gradient is 
more shallow and we do not see this effect. 

This demonstrates that, as expected, a gauge choice 
that leads to numerically stable evolutions with one evo- 
lution system may well be unstable for others. 



V. DISCUSSION 

We have investigated the numerical properties of the 
previously untested AA evolution system by a variety 
of numerical experiments, which for comparison we also 
performed for the ADM and BSSN systems. Since one 
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aspect of the test suite put forward in |2l| is that different 
numerical implementations of one and the same evolution 
system can be compared, let us point out that our results 
for ADM and BSSN a gree with |2ll for the specific resuhs 
shown there (see also [3l|). 

However, if the Courant factor in the robust stabil- 
ity test is lowered from 0.5 (as suggested in to 
0.25 we find that the stability of ADM dramatically im- 
proves (compare Figs. Q] and [21), while the stability prop- 
erties of AA remain unchanged. This implies that the 
Courant factor used in 21] is too large for the ADM sys- 
tem causing immediate exponential growth, instead of 
the initially linear growth expected for a weakly hyper- 
bolic system. This confirms a similar observation already 
described in . We use iterative Crank-Nicholson in all 
our simulations, but [24j also point out that dissipation 
can mask the linear growth expected for ADM. So one 
should repeat these experiments with a less dissipative 
scheme like third order Runge-Kutta to look for linear 
growth in ADM, but also in the AA systems, see Sec. Ill Al 

One property of BSSN that has not been explic- 
itly stated in |21| is its drastic failure on periodic do- 
mains. This observation is consistent with, and actually 
strengthens, the notion that BSSN performs well because 
it is able to propagate modes off the grid in partic- 
ular for isolated systems with radiative boundary condi- 
tions. It will be interesting to see how AA, ADM, and 
BSSN perform on a gauge wave test with outer boundary 
conditions that let the gauge wave propagate away from 
the domain. 

Our findings for the AA system can be summarized 
as follows. In the robust stability test, the AA system 
is stable for a long time, but eventually does go unsta- 
ble. In this case the A A system does not do as well as 
the BSSN system. For gauge waves in periodic geometry, 
the AA system is unstable, but the runs last much longer 
than the corresponding BSSN runs. The linear wave test 
shows that the AA system produces a larger drift in the 
phase than the ADM and BSSN systems, causing the 
linear wave to propagate faster compared to the analyt- 
ical solution in the AA evolution than in the other two 
evolutions. 

In the black hole runs, AA does as well as ADM and 



BSSN for geodesic slicing runs, but fails when trying 
to use gauges that makes BSSN runs long term stable. 
There are two distinct aspects of this gauge choice. While 
the gauge is appropriate geometrically (singularity avoid- 
ing and countering slice-stretching) , there is no reason to 
expect numerical stability for the AA system. Since it 
was non-trivial to find a numerically stable gauge choice 
for BSSN, it is not too surprising that additional work is 
required to find a gauge choice for the AA system that 
allows long run times for a single black hole, if indeed 
such a choice exists. 

We have found that the tests published in ^ are help- 
ful in exploring the properties of the AA system, but also 
that, not unexpectedly, these tests cannot by themselves 
determine whether an evolution system is worth explor- 
ing further for a particular physical system, say black hole 
evolutions. Note that one should expect that for differ- 
ent physical initial data sets different evolution systems 
are optimal, see e.g. ,32]. Since the relationship between 
choice of evolution system, gauge choice, outer boundary 
conditions, and the physical properties of the problem 
we are trying to simulate is complicated, it is not clear 
that sufficient understanding of the numerics needed to 
evolve binary black holes can be gained by singling out 
for example the evolution system and ignoring the other 
issues involved. Ultimately, if one wants to evolve black 
holes, one should evolve black holes. 

While the AA system has had some success in our nu- 
merical experiments, tests with black holes, proper outer 
boundaries, and appropriate gauge choices are needed to 
determine if the AA system has a future in numerical 
relativity. 
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FIG. 1: Robust stability test for the ADM system. The legend is: solid: p = 1, dotted: p — 2, dash-dot: p = 4, Loo (ham) 
is shown both as a function of time (left) and iterations (right). Note the almost perfect alignment when plotting versus the 
number iterations, which is the signature of a grid mode. 




FIG. 2: Robust stability test for the ADM system, with Courant factor 0.25. The legend is: solid: p — 1, dotted: p = 2, 
dash-dot: p = 4. We show Loo(ham) on the left and Loo{trK) on the right. The exponentially growing mode seen in Fig.^is 
not present. 
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FIG. 3: Robust stability test for the BSSN system. The legend is: solid: p = 1, dotted: p = 2, dash-dot: p = 4. 
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FIG. 4: Robust stability for the AA system, run until 1000 crossing times. The legend is: solid: p = 1, dotted: p = 2, dash-dot: 
p = 4. 
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FIG. 5: Robust stability for the AA system, run until 10000 crossing times. The legend is: solid: p = 1, dotted: p — 2, 
dash-dot: p = 4. We see that at late times the constraint violations are growing, which will probably cause a crash if we evolve 
long enough. 




FIG. 6: Gauge stability test for the AA system with a wave amplitude A = 0.1. The legend is: solid: p = 1, dotted: p = 2, 
dash-dot: p = 4. As expected, the metric converges at early times, but there is exponential growth and the run crashes after 
about 100 crossing times. 
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FIG. 7: Gauge stability test for the AA system with a wave amphtude A = 0.01. The legend is: solid: p = 1, dotted: p = 2, 
dash-dot: p = 4. We see that with this smaller amplitude the runs last 10 times longer but still crash eventually. 




FIG. 8: Gauge stability test for the BSSN system with a wave amplitude A = 0.01. The legend is: solid: p = 1, dotted: p = 2, 
dash-dot: p = 4. Both the metric and the extrinsic curvature tensor become non-convergent in a short time. 
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FIG. 9: Linear wave stability test for (from left to right) the AA system, the BSSN system, and the ADM system. We compare 
the numerical wave at 100 crossing times to the analytical solution. The legend is: solid: analytical solution, dotted: p = 1, 
dash-dot: p = 2, dash: p = 4. 
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FIG. 10: Linear wave stability test for (from left to right) the AA system, the BSSN system, and the ADM system. We show 
the L2 norm of the difference between the analytical and numerical values at different crossing times. The legend is: solid: 
p = 1, dotted: p = 2, dash-dot: p = 4. 
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FIG. 11: Geodesic slicing of Schwarzschild. Shown is the order of convergence cr of the conforms! metric component gzz for the 
AA system in the computational domain near the black hole for t = and t = 3.0M. We observe second order convergence 
near the black hole. 
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FIG. 12: Geodesic slicing of Schwarzschild. The first panel shows the analytical conformal metric component g^z as a function 
of z in the inner region of the domain (same region as in Fig. llljl . The analytical solution is shown at 3 different times: t = 1.55 
(solid), t — 3.0 (dotted) and t = 3.125 (dash-dot). The other 3 panels show the absolute value of the relative difference of the 
numerical and analytical solutions at the 3 different times. In these plots the AA solution is solid, ADM is dotted, and BSSN 
is dash-dot. 



